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Abstract: SMC^ (Chopin et al., 2013) is an efficient algorithm for sequential estimation and 
state inference of state-space models. It generates Ng parameter particles 0™, and, for each 
9”^, it runs a particle filter of size (i.e. at each time step, particles are generated in the 
state space A). We discuss how to automatically calibrate in the course of the algorithm. 
Our approach relies on conditional Sequential Monte Carlo updates, monitoring the state of the 
pseudo random number generator and on an estimator of the variance of the unbiased estimate 
of the likelihood that is produced by the particle filters, which is obtained using nonparametric 
regression techniques. We observe that our approach is both less CPU intensive and with smaller 
Monte Carlo errors than the initial version of SMC^. 
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1. INTRODUCTION 


Consider a state-space model, with parameter 9 G 0, la¬ 
tent Markov process (xt)t>oj and observed process {yt)t> 0 j 
taking values respectively in X and y. The model is 
dehned through the following probability densities: 9 has 
prior p{9), {xt)t>o has initial law pLg{xf) and Markov 
transition f^{xt\xt-i), and the pfs are conditionally inde¬ 
pendent, given the xfs, with density fj{yt\xt). Sequential 
analysis of such a model amounts to computing recursively 
(in t) the posterior distributions 


p{Q,xo,t\yo:t) = 
p{9)p,g{xo) 
p{y 0 :t) 


Y[fi^{Xs\Xs-l) 


n ivs 


.s=0 


or some of its marginals (e.g. p{9\yo:t))', the normalising 
constant p{yo:t) of the above density is the marginal 
likelihood (evidence) of the data observed up to time t. 


For a hxed 9, the standard approach to sequential analysis 
of state-space models is particle filtering: one propagates 
Nx particles in X over time through mutation steps (based 
on proposal distribution qtp{xt\xt-i) at time t) and resam¬ 
pling steps; see Algorithm I. Note the conventions: 1 : 
denotes the set of integers {1,..., yo:t is {yo,...,yt), 
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The output of Algorithm 1 may be used in different ways: 
at time t, the quantity J2n=i is a consistent 

(as Nx —>■ -boo) estimator of the filtering expectation 
¥\Lp{xt)\yci..t,9\, In addition, tt(d) is an unbiased estimator 
of incremental likelihood p{yt\yo:t-i,d), and ns=oh^) 
is an unbiased estimator of the full likelihood p{yo:tW) 
(Del Moral, 1996, Lemma 3). 

In order to perform joint inference on parameter 9 
and state variables, Chopin et al. (2013) derived the 
SMC^ sampler, that is, a SMC (Sequential Monte Carlo) 
algorithm in 0—space, which generates and propagates Ng 
values 6*™ in 0, and which, for each 6*™, runs a particle 
filter (i.e. Algorithm 1) for 9 = 9™, of size N^. One issue 
however is how to choose N^: if too big, then CPU time is 
wasted, while if taken too small, then the performance of 
the algorithm deteriorates. Chopin et al. (2013) give formal 
results (adapted from Andrieu et al. (2010)) that suggest 
that Nx should grow at a linear rate during the course of 
the algorithm. They also propose a practical method for in¬ 
creasing Nx adaptively, based on an importance sampling 
step where the Ng particle systems, of size N^, are replaced 
by new particle systems of size Nf™. But this importance 
sampling step increases the degeneracy of the weights, 
which in return may leads to more frequent resampling 
steps, which are expensive. In this paper, we derive an 



Algorithm 1. Particle filter (PF, for fixed 8 ) 


Algorithm 2. IBIS 


Operations involving superscript n must be performed for 
all n G 1 : N^. 

At time 0: 

(a) Sample Xq ^ qo,9ixo). 

(b) Compute weights 

... ^^e{x1^)f^{yo\xo) 

U’ 0 ,e(XQ) — „(rr.n\ 

qofiyxQ) 

normalised weights, W^g = wo,e(a;o)/'*«o,e(aio), 
and incremental likelihood estimate 

U0) = N-^Y.n=i<,e- 

Recursively, from time < = 1 to time t = T: 

(a) Sample a” ~ the multinomial distribu¬ 
tion which generates value i G 1 : with probability 

WUe- 

(b) Sample xj* ~ 

(c) Compute weights 


wt,e{xt^i,xt) = 




and incremental likelihood estimate 

m = NAj:nAm,9i<A,A)- 


Operations involving superscript m must be performed for 
all m G 1 : N 9 . 

(Init) Sample d'" ~ set w™ ^ 1. 

From time t = 0 to time t = T, do 
(a) Update importance weights 

w™ ^ w™ X p{yt\yo-.t-i,d)- 


(b) If ESS< ESSmin, sample (for all m) 0™ from 
mixture 

, Ne 

^ 75 — 


V ^ 


m—1 


where Kt is a Markov kernel with invariant distribution 
p(d|?/ 0 :t); finally reset particle system to 

Ql-.Ne ^ gV.Ne ^ ^ (1, . . . , 1). 


to the spread of the current particle population (i.e. vari¬ 
ance of random step equals some fraction of the covariance 
matrix of the resampled particles). 

The main limitation of IBIS is that it requires evaluating 
the likelihood increment p{yt\yQ:t-i,Q)) which is typically 
intractable for state-space models. On the other hand, we 
have seen that this quantity may be estimated unbiasedly 
by particle hltering. This suggests combining IBIS (i.e. 
SMC in the d-dimension) with particle filtering (i.e. SMC 
in the x*—dimension), as done in the SMC^ algorithm. 


alternative way to increase adaptively, which is not 
based on importance sampling, but rather on a CSMC 
(conditional Sequential Monte Carlo) update, which is less 
CPU intensive. 

2. BACKGROUND ON SMC^ 


2.1 IBIS 


To explain SMC^, we first recall the structure of the IBIS 
algorithm (Chopin, 2002) as Algorithm 2. For a model with 
parameter 0 G 0, prior p{0), data t/o:T: and incremental 
likelihood p{yt\yo:t-i,0)^ IBIS provides at each iteration t 
an approximation of partial posterior p(d|?/o:t)- In practice, 
IBIS samples Ng particles 0™ from the prior, then perfoms 
sequential importance sampling steps, from p(d|j/o:i-i) to 
p{0\yo-.t) using incremental weight p{0\yt),t)/p{ 0 \yt 3 -t-i) oc 

p{yt\vo:t-iA)- 


To avoid weight degeneracy, one performs a resample-move 
step (described as Step (b) in Algorithm 2). When the ESS 
(effective sample size) of the weights, computed as: 


(Sr-Ne rn\2 

ESS(wi^^«) = - - G [1,A^] 

goes below some threshold ESSmin (e.g. N/ 2 ), the d^’s 
are resampled, then moved according to some Markov 
kernel Kt that leaves invariant the current target of the al¬ 
gorithm, p{9\yti.,t). This resample-move step re-introduces 
diversity among the 0-particles. 


A convenient default choice for Kt is several iterations of 
random-walk Metropolis, with the random step calibrated 


2.2 SMCP 


The general structure of SMC^ is recalled as Algorithm 3. 
Essentially, one recognises the IBIS algorithm, where the 
intractable incremental weight p{yt\yo:t-i,0"^) has been 
replaced by the unbiased estimate £((0™). This estimate is 
obtained from a PF run for 0 = 0™; thus N 9 PFs are run 
in parallel. Denote the random variables 

generated by the PF associated to 0™. 

This ‘double-layer’ structure suggests that SMC^ suffers 
from two levels of approximation, and as such that it 
requires both —>• -too and N 9 —)• -foo to converge. 
It turns out however that SMC^ is valid for any fixed 
value of Nx] that is, for any fixed > 1, it converges as 
N 9 —>■ -foo. 


This property is intuitive in the simplified case when 
resampling-move steps are never triggered (i.e. take 
ESSmin = 0)- Then SMC^ collapses to importance sam¬ 
pling, with weights replaced by unbiased estimates, and it 
is easy to show convergence from first principles. 


We now give a brief outline of the formal justification of 
SMC^ for fixed N^, and refer to Chopin et al. (2013) for 
more details. SMC^ may be formalised as a SMC sampler 
for the sequence of extended distributions: 
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where ipt,e denotes the joint pdf of the random variables 
generated by a PF up to time t (for parameter 0), and 
C(0) denotes the unbiased estimate of the likelihood incre- 



Algorithm 3. SMC^ 


Operations involving superscript m must be performed for 
all m € 1 : 

(Init) Sample 0”^ ~ p(^)) set w™ ^ 1. 

From time t = 0 to time t = T, do 

(a) For each fl'", run iteration t of Algorithm 1, so as to 

obtain and 

(b) Update weights 

(c) If ESS(a;^'^®) < ESSmin, sample (for all m) 

from mixture 


1 


L^rn —1 


Ne 

■E 

m—l 




l-.Nxim l:Nx,m 


l:t 


^\d- 
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where Kt is a PMCMC kernel with invariant distribu¬ 
tion Trt{9,XQ.^'^ 
system to 
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ment computed from that PF, ioi9) = J2n=i ^oi^^), 
4(6») = N-'^J 2 n=iWs, 0 ix'',ti,x^) for s > 0; i.e. 4(6») is 
actually a function of {9,Xq.^'‘ 

One recognises in tt^ the type of extended target distri¬ 
bution simulated by PMCMC (Particle MCMC, Andrieu 
et al. (2010)) algorithms. Note irt is a proper probabil¬ 
ity density (it integrates to one), and that the marginal 
distribution of 9 is p(0|j/o:i)- These two properties are 
easily deduced from the unbiasedness of n!=o^s(^) 
an estimator oi p{yo-,t\9)). In addition, 


7r,(0,4E,a;:E) = 


l:Afx 


7rt_i(6»,xJ:^"i,ai:(E) 
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where one recognises in the second factor the distribution 
of the variables generated by a PF at time t, conditional 
on those variables generated up to time t — 1. Thus, the 
equation above justifies both Step (a) of Algorithm 3, 
where the particle filters are extended from time t — 1 to 
t, and Step (b), where the particles (0™, 
are reweighted by ^((0™). 


We describe in the following section PMCMC moves 
that may be used in Step (c). Before, we note that a 
naive implementation of SMC^ has a O{tNxN 0 ) mem¬ 
ory cost at time t, as one must stores in memory 
(0™, for each to G 1 : Ng. This memory 

cost may be substantial even on a modern computer. 


2.3 PMCMC moves 

To make more explicit the dependence of the unbiased es¬ 
timate of the likelihood on the variables generated during 
the course of PF, define 


Tt(6',4:rAonf")= 

'' Nx t r Nx ^ ') 

E ^o.eixo) f n 1 AT E '^s,eix‘'/_i,x^) I . 


1 

iW 


The PMMH (Particle Markov Metropolis-Hastings) ker¬ 
nel, described as Algorithm 4, may be described infor¬ 
mally as a Metropolis step in 0-space, where the likeli¬ 
hood of both the current value and the proposed value 
have been replaced by unbiased estimators. Formally, as 
proven in Andrieu et al. (2010), it is in fact a standard 
Metropolis step with respect to the extended distribu¬ 
tion 7rt(0, in particular it leaves invariant 

pi9\yo:t)- (For convenience, our description of PMMH as¬ 
sumes a random walk proposal, but PMMH is not re¬ 
stricted to this kind of proposal.) 


Algorithm 4. Random walk PMMH update 


Input: (0,4:E,ai;E) 

Output: (0,iJ:fAhi[:E) 

1. 9* = 9 +z, lV(0,St). 

2. Generate PF (Algorithm 1) for parameter 0*; let 

the output. 

3. With probability 1 A r. 


p(0*)L,(0Aa:i:E’*,ai:E”^) 


l-.Nx 

o-.t 

l:Nx 

0:t 


p(0)L,(0,xi:E,a}:^) 


let (0,iJ;E,ai:E) 
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) (0, Xq] 
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l-.N, 

t 


.ai:t 




otherwise 


In practice, we set the covariance matrix of the 
proposal, to a fraction of the covariance matrix of the 
resampled 0-particles. 

One advantage of using PMHMH within SMC^ is that 
it does not require storing all the variables generated 
by the Ng PFs: operations at time t > 0 require only 
having access to, for each to, (0™, 0 ( 1 ^“'’™) and 

Lt_i(0™, which is computed recursively. 

Memory cost then reduces to 0{NgNx). 

The Particle Gibbs approach is an alternative PMCMC 
step, based on the following property of target tt*: if one 
extends ttj with random index fc, such that k € 1 : N^, 
and k ~ A4(1F^'^“'), the normalised weighs at the final 
iteration, then (a) the selected trajectory, together with 
0, follow the posterior distribution p(0, a;o:t|j/o:t)i (b) 
the remaining arguments of tt* follow a CSMC (conditional 
SMC) distribution, which corresponds to the distribution 
of the random variables generated by a PF, but conditional 
on one trajectory fixed to the selected trajectory; see 
Algorithm 5. 

In contrast with PMMH, implementing particle Gibbs 
steps within SMC^ requires having access to all the vari¬ 
ables (0™, at time t, which as we have 

already discussed, might incur too big a memory cost. 



Algorithm 5. Particle Gibbs update 

Input: 

Output: , 0 , 1 '.^'°) 

1. Sample h ~ with W*" 

From s = < - 1 to 

s = 0, set hs ^ Og+i • Set x], ^ x^g% a]. = 1 for all 
s€0 :T.^ 

2 . Sample 9 from a MCMC step that leaves invariant 
distribution p(0|a:o:t, 2 /o:t)) but with xo-,t set to Xq.^. 

3. Sample as in Algorithm 1, but for pa¬ 

rameter 9 and conditionally on Xq.^., that is: at time 0, 
generate ig ~ ^ for n € 2 : A^, at time 1, sample 

a” ~ for n € 2 : iV, and x” ~ 

and so on. 


2 . 4 ^ Choosing 

Andrieu et al. (2010) show that, in order to obtain reason¬ 
able performance for PMMH, one should take = 0{t). 
Andrieu et al. (2013) show a similar result for Particle 
Gibbs. 

In the context of SMC^, this suggests that should 
be allowed to increase in the course of the algorithm. 
To that effect, Chopin et al. (2013) devised an exchange 
step, which consists in exchanging the current particle 
systems, of size N^, with new particle systems, of size 
Al)?®'", through importance sampling. In Chopin et al. 
(2013)’s implementation, the exchange step is triggered 
each time the acceptance rate of the PMMH step (as 
performed in Step 3. of Algorithm 4) is below a certain 
threshold, and Al)?®" = 2Alr (i.e. doubles every time). 

The main drawback of this approach is that it introduces 
some weight degeneracy immediately after the resampling 
step. In particular, we will observe in our simulations that 
this prevents us from changing too frequently, as the 
ESS of the weights then becomes too low. 

In this paper, we discuss how to use a Particle Gibbs step 
in order to increase without changing the weights. 

3. PROPOSED APPROACH 

3.1 Particle Gibbs and memory cost 

'We first remark that the Particle Gibbs step. Algorithm 5, 
offers a very simple way to change during the course of 
the algorithm: In Step (2), simply re-generate a particle 
system (conditional on selected trajectory ig.J of size 
A^"®". But, as already discussed, such a strategy requires 
then to access past particle values a;" (and also a"), rather 
than only current particle values x". 

This problem may be addressed in two ways. First, one 
may remark that, to implement Particle Gibbs, one needs 
to store only those x" (and a”) which have descendant 
among the current particles x”. Jacob et al. (2013) 
developed such a path storage approach, and gave con¬ 
ditions on the mixing of Markov chain (x*) under which 
this approach has memory cost 0{t + N^logN^) (for a 


single PF with particles, run until time t). Thus, an 
implementation of this approach within SMG^ would lead 
to a 0{Ng{t + Nx log A^j;)) memory cost. 

A second approach, developed here, exploits the determin¬ 
istic nature of PRNGs (pseudo-random number genera¬ 
tors): a sequence zg, zi,..., Zj,... of computer-generated 
random variates is actually a deterministic sequence de¬ 
termined by the initial state (seed) of the PRNG. It is 
sufficient to store that initial state and Zg in order to 
recover any Zi in the future. The trade-off is an increase 
in CPU cost, as each access to Zi require re-computing 

Zi, . . . ,Zi. 

We apply this idea to the variables By 

close inspection of Algorithm 3, we note that variables in 
a ‘time slice’ (xi’^-™, 0 < s < t (or Xg^"’™ at 

time 0) are always generated jointly, either during Step 
(a), or during Step (c). In both cases, this time-slice is 
a deterministic function of the current PRNG state and 
the previous time slice. Thus, one may recover any time 
slice (when needed) by storing only (i) the PNRG state 
(immediately before the generation of the time slice); and 
(ii) in which Step (either (a) or (c)) the time slice was 
generated. This reduces the memory cost of SMC^ from 
OitNgN,) to 0{Ne{t + N,)). 

Compared to the path storage approach mentioned above, 
our PRNG recycling approach has a larger CPU cost, a 
smaller memory cost, and does not require any conditions 
on the mixing properties of process (x(). Note that the 
CPU cost increase is within a factor of two, because each 
time a Particle Gibbs update is performed, the number 
of random variables that must be re-generated (i.e. the 
x" and a" in Algorithm 5) roughly equals the number of 
random variables that are generated for the first time (i.e. 
the x” and a” in Algorithm 5). 

3.2 Nonparametric estimation of 

As seen in Algorithm 3, a Particle Gibbs step will be 
performed each time the ESS goes below some threshold. 
That the ESS is low may indicate that is also too low, 
and therefore that the variance of the likelihood estimates 
Lt(0™, is too high. Our strategy is to 

update (each time a Particle Gibbs step is performed) 
the current value of to A"®" = r/d^, where is 
some (possibly rough) estimate of the variance of the 
log likelihood estimates. This is motivated by results 
from Doucet et al. (2012), who also develop some theory 
that supports choosing r « 1 is optimal (although their 
optimality results do not extend straightforwardly to our 
settings). 

Assume 0 C M”*. To estimate we use backfitting to 
fit a GAM (generalized additive model) to the responses 
K - log Lt[0 ). 

d 

R'" = a + ^/,(G7)+£”^, 

i=i 

using as covariates (7™ the d principal components of 
the resampled 6*-particles. The estimate is then the 
empirical variance of the residuals. See e.g. Chap. 9 of 
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Fig. 1. Marginal posterior p(cr^, p|j/o:i 5 )) as approximated 
by SMC^ run until t = 15, and linearly transformed 
so that axes are the two principal components. 

Hastie et al. (2009) for more details on backfitting and 
GAM modelling. 

We found this strategy to work well, with the caveat that 
choosing r required some trial and error. 

3.3 Additional considerations 

Using Particle Gibbs as our PMGMG move within 
SMG^ hast two advantages: (a) it makes it possible to 
change without changing the weights, as explained 
above; and (b) it also makes it possible to update the 
0"* according to Gibbs or Metropolis step that leaves 
d\xo:t,yo:t invariant); see Step (3) of Algorithm 5. For 
models where sampling from 9\xo-,t,yo:t is not convenient, 
one may instead update 9 through several PMMH steps 
performed after the Particle Gibbs step. 

4. NUMERICAL EXAMPLE 

We consider the following stochastic volatility model: xq ^ 
N(/r,(T2/(l-p2)), xt-y = p{xt-i-y.)+aet, et N{0, 1) 
and yt\xt ~ N{0,e^^)] thus 9 = {p,p,a), with p € [—1,1], 
cr > 0. We assign independent priors to the components 
oi 9: p ^ N(0, 2^), p ^ N{0, 1) constrained to [—1,1], and 
~ /G(3, 0.5). The dataset consists in log-returns from 
the monthly SP500 index, observed from 29/05/2013 to 
19/12/2014; T = 401. 

Figure 1 plots the marginal posterior p(p, cr^lyons), as 
approximated by SMC^, run up to time 15. This figure 
illustrates the need for modelling nonparametrically the 
true likelihood as a function of 9, in order to estimate the 
variance of the estimated likelihood. 

For this model, sampling jointly from 9\xo.t,yo:t is difficult, 
but it is easy to perform a Gibbs step that leaves invariant 
9\xo-.t,yo-.t, as the full conditionals of each component (e.g. 


Fig. 2. Evolution of over time for 5 runs of the 
four considered SMC^ algorithms; red dotted line is 
Algorithm (a), blue dashed is (b), black solid is (c), 
green double-dashed is (d). Results of (c) and (d) are 
nearly undistinguishable. 

p\a, p, xo:t, yo-.t and so on) are standard distributions. Let’s 
call ‘full PG’ Algorithm 5, where Step 2 consists of this 
Gibbs step for 9\xo.t, yo-.t', and conversely let’s call ‘partial 
PG’ Algorithm 5 with 9 = 9 in Step 2 {9 is not updated). 

We compare four versions of SMC^: (a) the standard 
version, as proposed in Chopin et al. (2013) (i.e. Step (c) 
of Algorithm 3 is a PMMH step, and that step is followed 
by an exchange step to double when the acceptance 
rate of PMMH is below 20%); (b) the same algorithm, 
except that an exchange step is systematically performed 
after Step (c), and is set to the value obtained with 
our non-parametric approach (see Section 3.2); (c) the 
version developed in this paper, with full PG steps (and 
Nx updated through the non-parametric procedure); (d) 
the same algorithm, but with partial PG steps, followed 
by 3 PMMH steps to update 9. 

The point of Algorithm (b) is to show that adapting too 
often during the course of the algorithm is not desirable 
when using the exchange step, as this leads to too much 
variance. The point of Algorithm (d) is to see how our 
approach performs when sampling from 9\xo,t, yo-.t (either 
independently or through MCMC) is not feasible. 

Figure 2 plots the evolution of over time for the four 
SMC^ algorithms. One sees that, for these model and 
dataset, the CPU cost of the standard SMC^ algorithm 
is quite volatile, as increases very quickly in certain 
runs. In fact certain runs are incomplete, as they were 
stopped when the CPU time exceeded 10 hours. On the 
other hand, the CPU cost of other versions is more stable 
across runs, and, more importantly, quite lower. 

Figure 3 plots the empirical variance of the estimated 
marginal likelihood (evidence, p{yo:t)), normalised with 
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Fig. 3. Empirical variance of estimated marginal likelihood 
piyo:t) multiplied by average CPU time; same legend 
as Figure 2, results from Algorithm (b) are omitted. 
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Fig. 4. PMMH acceptance rate across time; same legend 
as Figure 2. Black line marks 20% target. 

the running time up to time step t. One observes that 
version (c) does quite better than (d), and far much better 
than (a). Results from Algorithm (b) were to variable to 
be included. 

Figure 4 plots the acceptance rate of PMMH steps for 
Algorithms (a), (b) and (d). (Recall that Algorithm (c) 
does not perform PMMH steps). Note the poor perfor¬ 
mance of Algorithm (b). Figure 5 compares the box-plots 
of posterior estimates of a at hnal time T, obtained from 
several runs of Algorithms (c) and (d). Algorithm (c) 
shows slightly less variability, while being 30% faster on 


Fig. 5. Box-plots of posterior estimate of parameter a at 
hnal time T, over repeated runs of Algorithm (c) (left 
panel) and Algorithm (d) (right panel). 

average. One sees that the improvement brought by ability 
to sample from 0 \xo-,t,yo-.t is modest here for parameter 
estimation, but recall that in Figure 3, the improvement 
was more substantial. 

ACKNOWLEDGEMENTS 

We thank Pierre Jacob for useful comments. 

REFERENCES 

Andrieu, C., Doucet, A., and Holenstein, R. (2010). 
Particle Markov chain Monte Carlo methods. J. R. 
Statist. Soc. B, 72(3), 269-342. doi:10.1111/j.l467- 

9868.2009.00736.x. 

Andrieu, C., Lee, A., and Vihola, M. (2013). Uniform Er- 
godicity of the Iterated Conditional SMC and Geometric 
Ergodicity of Particle Gibbs samplers. ArXiv e-prints. 
Ghopin, N. (2002). A sequential particle hlter for static 
models. Biometrika, 89, 539-552. 

Ghopin, N., Jacob, P., and Papaspiliopoulos, O. (2013). 
SMG^: A sequential Monte Garlo algorithm with parti¬ 
cle Markov chain Monte Carlo updates. J. R. Statist. 
Soc. B, 75(3), 397-426. 

Del Moral, P. (1996). Non-linear hltering: interacting 
particle resolution. Markov processes and related fields, 
2(4), 555-581. 

Doucet, A., Pitt, M., Deligiannidis, G., and Kohn, R. 
(2012). Efficient implementation of Markov chain Monte 
Carlo when using an unbiased likelihood estimator. 
ArXiv preprint. 

Hastie, T., Tibshirani, R., Friedman, J., Hastie, T., Fried¬ 
man, J., and Tibshirani, R. (2009). The elements of 
statistical learning, volume 2. Springer. 

Jacob, P., Murray, L., and Rubenthaler, S. (2013). 
Path storage in the particle hlter. Statist. Corn- 
put, 1-10. doi:10.1007/sll222-013-9445-x. URL 

http;//dx.doi.org/10.1007/s11222-013-9445-x. 























